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ABSTRACT 



The effects of trim on stability of motion during depth 
control of submersible vehicles are analysed. Full state 
feedback control is used to provide stable response in the dive 
plane, and feedforward control is used to ensure steady state 
accuracy. A complete set of stability maps is generated for 
various values of metacentric height, longititunal center of 
gravity/center of buoyancy separation, forward speed, and 
control law time constant. The results clearly indicate ranges 
of parameters that should be chosen in design and operation of 
a given vehicle. 
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NOMENCLATURE 



Symbol : 
A 



a 

b(x) 

B 

B 

C D 

C l r c 2 

d q '<K 

5 b 

6 S ,6 

Inrlp 

k 1 ,k 2 ,k 3/ k 4 

K 0 

m 

M 

M a 

0 

q 

t c 

u 

Uo 

w 

w 

x 

^B / Z B 



r Z G 
Z 




Definition ; 

Closed loop dynamics matrix for 
linearized system 

Control surface coordination gain 

Local beam of the hull 

Vehicle buoyancy 

Control matrix in state space 

Quadratic drag coefficient 

Coupled heave and pitch terms 

Cross flow drag terms 

Bow plane deflection 

Stern plane deflection 

Vehicle mass moment of inertia 

Cross flow drag terms 

Controller gains in 9,w,q, and z 

Feedforward gain 

Vehicle mass 

Pitch moment 

Partial derivative of M w.r.t. a 

Pitch angle 

Pitch rate 

Time constant 

Forward speed 

Nominal speed 

Heave velocity 

Vehicle weight 

State vector 

Body fixed coordinates of vehicle 
center of buoyancy 

Body fixed coordinates of vehicle 
center of gravity 
Deviation off the nominal depth 
Vehicle metacentric height 
Heave force 

Partial derivative of Z w.r.t. a 



the 



vi 



I . INTRODUCTION 

The fundamental goal of submarine control is to reach and 
maintain ordered depth. Experimental designs involve expensive 
model testing such as Darpa Suboff Model (DTRC Model 5470) 
[Ref. 6]. Much research has been done in depth control of 
submarines [Ref. 3,5]. Our goal is to develop an analytic 
method to determine the stability properties of a design. 

The stability of a design will have a significant impact 
on its responsiveness. A vehicle with a large margin to 
instability will not be very responsive. The problem becomes 
one of determining how close to the margins we can get without 
a total loss of stability. By expanding the scope of our 
research to include nonlinear terms we are able to define the 
limits of stability and therefore margins. 

Previous studies analyzed stability properties of the 
system, specially static bifurcations [Ref. 2] and 
bifurcations to periodic solutions [Ref. 1]. The latter study 
which is used as a basis for this work, was restricted to 
level, zero trim flight paths. 

The purpose of this thesis is to develop a program for 
finding the limits of stability for an out of trim submarine 
at moderate and high speeds. These limits are mapped using a 
Hopf bifurcation analysis program included in the Appendix. 
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II. EQUATIONS OF MOTION 



The motion of the submersible in the vertical plane can be 
modeled by four coupled nonlinear differential equations for 
pitch rate (q) , heave velocity (w) , pitch angle (0) and heave 
( z ) . With a body fixed coordinate frame at the vehicle * s 
geometric center , we can express Newton's equations of motion 
as 




( 2 . 1 ) 




l 



(w-xq) 3 



I y - mx G (w - Uq ) -z G wqm = -x gb Bcos 0 - z GB Bsin 0 b s 




(2.3) 



z - - Usin 0 -mv cos 0 



(2.4) 
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Equations 2 • 1 through 2 . 2 can be written in a more compact 



form as, 



w = a. 



,Uw 



M a 






b 2 U 2 b b + d w (w,q) + c 1 (w.q) 



(2.5) 



q = a 21 Uw + a 22 Uq + a^z^sin 0 + a 23 x QB cos 0 + bjU 2 t> s 
b 2 U 2 b h + d (w,q) + c 2 (w,q) 



( 2 . 6 ) 



where, 

D v = ( m ~ ( mx G + Z 4>( mX G + M J 

djjD v = (Z„-2C D E 0 Utan Q g ) (I y - M ^) + (mx G + Z (M 2 C ^E 2 U tan Qq) 

a i2 D v = ( m+z q + 2C D EjUtan Q 0 )(I y -M^) + 

(mx G + Zq) (M q - mx Q - mz Q U tan d Q - 2 C D E 2 U tan d Q ) 

a 21 D v = (M w + 2C D E I Utan d Q )(m - ZJ + (mx G + Mty(Z w ~ 2C D E Q Utan Q g ) 

a 22 D v = (M q - mx q - mz G U tan 0 O - 2 C D E 2 U tan 0 Q ) 

(mx G + M^) (m + Z q + 2 C D E X t/tan 0 O ) 

(2.7) 

a 23 D v = ~( m + Z JB 
a 13 D v = " ( mX G + 

b iD v - (I y -M 4 )Z,-( G +Z 4 )M h 
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b 2 D v = ( m ■ Z J M b + (mx G + MJ 



d q (w.q)D v = (m - Z^/ ? + (mx G + M*)I w 



dJw,q)D v = (I y -M 4 )I„- (mx G + 1J1 q 



Cj(w,q)D v = (I y - M 4 )m z G q 2 - (mx G + ZJmz^q 



c 2 (w.q)D v = - (m - Z ylr )mz G wq ~ (mx G + M J m z Q q 2 



In these equations the submersible is assumed to be 
neutrally buoyant (W=B) , and statically stable ( z G > z B ) . Here 
we can assume z B to be zero, hence z GB = z G 

At steady state the cross flow drag integral terms I H and 
I p have the form, 



From equation (2.3) it is seen that w is equal to tan0 o at 
steady state. The fb(x)dx term is computed numerically for the 
SUBOFF model as E 0 , and Jb(x)xdx term as E u Therefore, the 
cross flow drag terms become, 





( 2 . 8 ) 



I H = -C B w /W/ E 0 



I P = C D w /w/ E 



(2.9) 
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Because we have two rudders at the bow and the stern, our 
system of equations is multi-input. To reduce this system into 
a single input system the linear combination of the control 
inputs will be modified into the following form, 

6 = 5 S , 5 b = ad s (2.10) 
where a is the control surface coordination gain. The value 
of a ranges from -1 to 1. The selection of the value of a will 
allow the planes to operate as desired for the particular 
maneuvering condition, i.e., a = 0 for no bow plane control, 
a = -1 for bow plane and stern plane opposed to each other, 
yielding the maximum pitch moment, and a = 1 for bow and stern 
plane control in the same direction, yielding the maximum 
heave force. 
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FIGURE 2.1 Geometric representations of the basic 



definitions . 
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CONTROL LAW 



A. INTRODUCTION 

The control design problem can be expressed in state space 
as follows, 

. , n* (3.1) 

x = Ax + B 6 ' ' 

where the state vector is 



x = 



6 

w 

q 



(3.2) 



Equation 3.1 in our case is. 



e 




0 


0 


1 


0 




0 






a 13 Z GB~ b l u2k l 


djjti - bjU k 2 


a jj4 - bjU k s 


-bp 2 k 4 




w 


q 




a 23 Z GB ~ b i l 2fc l 


a 2} u -bji 2 k 2 


a 2 ji ~ b ? 2 k 3 


-bji 2 k 4 




q 


Z j 




-u 


1 


0 


0 




z 



Our aim is to find a controller which will assure us a 
stable closed loop system. The only control input is the dive 
plane angle, 5. 
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B . FEEDBACK CONTROL 



1 . Pole Placement 

The full state feedback controller is a linear function of 
the states and has the form, 

5 = - K x (3.4) 

where K is the vector of feedback gains which are to be 
determined in order to give the desired closed loop system 
dynamics. Substituting equation 2.12 into 2.10 yields, 

x = (A- BK )x (3.5) 

The feedback gains K must be chosen such that A - BK has the 

desired eigenvalues . The actual characteristic equation of the 
closed loop system is given by, 

det (A - BK - si) = 0 (3.6) 

The required values of K are obtained by matching coefficients 
in the two polynomials of the actual and the desired 
characteristic equations. Equation 3.5 becomes. 



e 




0 


0 


1 


0 




0 




0 


vi» 


__ 


a 13?GB 


a n u 


a l? 


0 




w 


+ 


bjU 2 


4 




a 2?GB 


a 2l“ 


a 2? 


0 




<1 




b 2 u 2 


z 




-u 


1 


0 


0 




z 




0 



The characteristic equation of the 



closed loop system is. 
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del 



-s 0 10 

2 2 2 2 
a lS z QQ-b J u kj a ]I u -bjU k 2 -s arfi-bjU ^3 ~bj u * 4 

a 2?GB~ b i l2k l a 21 U ~ b ^ 2fc 2 a 2^~ b ^ 2k 3~ s ~ b 4 

l 0 -s 



-u 



(3.8) 



which reduces to. 



where, 



( A 2 k 2 + 


A s k - 


i-Ej) 


° 3 u- 


-Bjkj B 2 k 2 E 3 k 3 b 4 k 4 


(-c 


Si- 


C 2 k 2 


-C 3 k 3 


- E 3 )s + (-D I k 4 - D 2 k 4 ) = 


= 


ll 2 








= b 2 


u 2 








b i ~ 


a 12 


b 2 ) 


u 3 




: ( a u 


b 2 


~ a 2i b i) 


u 3 


’■ ( a 23 


b> 


~ a i3 


b 2 ) 


z gb U 2 


b ! + 


b 2 


~ a i2 b 2 ) 


u 3 


b 2 " 


a 21 


b 1 ) 


u 4 





(3.9) 



(3.10) 



Ei = (a n + a 22 ) u 

E 2 = a 2J z GB + (a 12 a 21 - a n a 22 ) u 2 

E 3 = ( a i3 a 21 ~ a il a 23 ) Z GB U 

Now, let's assume that we want to place the closed loop poles 
at -pj, -p 2 , -p 3 , -p 4 to have the desired system response. Then 
the desired characteristic equation is, 

432 n ( 3 * 11 ) 

s + <tjS + (t 2 s + tt^s + a. 4 = 0 ' ' 
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where , 



a i = Pi + P 2 + P 3 + P 4 

= Pi P 2 + Pi P 3 +Pi P< + P2P3 + P2P4 + P3P4 (3.12) 

<*3 = Pi P2 P3 + Pi P2 P 4 + Pi P3 P 4 + P2 P3 P 4 
<** = Pi P2 P 3 P 4 

The feedback gains can now be computed by equating the 
coefficients of equation 3.9 and 3.11, 

A 2 k 2 + A 3 k 3 = -a 2 - E 1 

B 1 k 2 + B 2 k 2 + B 3 k 3 + B 4 k 4 = a 2 + E 2 (3.13) 

C 3 k 3 + C 2 k 2 + C 4 k 4 = a 3 + E 3 
(D, + D 2 ) k 4 = cr 4 

We established the method for placing the poles of the 
system, but we also need to know the desired locations of the 
poles . 

2 . Pole Location Selection 

In a typical second order system control law design , 
transient response specifications are given. This results in 
an allowable region in the s-plane where the desired location 
of the poles can be obtained. For higher order systems the 
concept of dominant roots can be employed. In selecting poles 
the physics of the system must be considered. If the poles are 
too negative, a small time constant will result, and the 
system may not be able to react that fast. If we use big gains 
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K, this also means that the control effort will be large. In 
practice there are limits based on actuator size and 
saturation . 

Considering the control law design to stabilize the 
submarine to a level flight path at 0 = 0 it is required that 
the submarine return to the level flight, after some small 
disturbances in 0 or z, within the time it takes for the 
vehicle to travel three ship lengths. Since our model is 14 
feet long and its velocity is 5 feet/sec., the required 
recovery time is about 10 seconds. This means the time 
constant is 3 seconds and the closed loop poles should be 
placed at approximately -0.3. 

After placing the poles using equation 3.13 the 
control law is found to be, 

5 = 0.9917 6 + 0.8333 w + 0.6026 q - 0.0351 z (3.14) 

C . FEEDFORWARD CONTROL 

The previous discussion on feedback controller assures 
closed loop stability, but it acts as a regulator in other 
words takes all the states to zero values. If we have constant 
disturbances or we want to track some reference value other 
than zero we can not do this with state feedback alone. 
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FIGURE 3.1 Feedback and feedforward control application to 
our linearized model 



In the case of non-zero set points or constant 
disturbances we again need to have the exact same full state 
feedback to ensure closed loop stability. But we also need to 
introduce an additional term to our controller in the form 

u = - K x + k 0 (3.15) 
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where k 0 is the constant feedforward term. k 0 is given as 
[Ref. 4], 

k 0 = H'^O) x 0 (3.16) 

where x 0 is the reference values of the states and H c _1 (s) is 
the closed loop transfer function, 

Hi'fs) = C (si - A + BK)' 1 B (3.17) 

Another way of getting k 0 is looking at the steady state 
equations of motion. In steady state all the time derivatives 
in equations 2.1 through 2.4 go to zero and we have. 



iv = tan 0 O 



Z b b+Z „ W ~ I H = 0 



(3.18) 

(3.19) 



~(x G - x b )Bcos Q q - z GB Bsin 0 o + A/ ft 8 +M w w + I p = 0 



(3.20) 



If the equation 3.19 is multiplied with M 6 and equation 3.20 
with Z 5 and set equal to each other and plug in the equation 
3.18, we have an equation depending only on 0. 

( Z y^U- M y/i) tan0 o + x g bz i cosQ o + z o BZ i sinQ o + 91 . 

C^EoM, - EjZ h ) tan Q 0 \tanQ 0 \ ( * ^ ^ ' 

Where E 0 and E x are the integral terms computed numerically for 
the Suboff model. The steady state value of 9 0 is 
found from equation 3.20 by using a Newton-Raphson method 
in Bifurl program in the Appendix. Then we can get 5 from 
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equation 3.19, 



6 = 




-C D E 0 tan Qj\tan 8 „| - Z w w 



(3.22) 



After getting 5, we easily find k 0 from equation 3.15, 

k 0 = <5 + Kx (3.23) 

A plot of the steady state angle 0 O versus x G for different 
values of z G is shown in Figure 3.2. 



14 




FIGURE 3.2 



Steady state pitch angle 0 O as metacentric 



height varies 
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IV. 



BIFURCATION ANALYSIS 



A. STABILITY 

The nonlinear equations of motion in the dive plane 2 . 1 
through 2.4 can be expressed in a compact form as follows, 



Where x is the state variable vector x=[0,w,q,z]. We know 
that the equilibrium points, x 0 of the system are defined by, 



This is a nonlinear system of algebraic equations and it 
may have multiple solutions in x 0 , which means that the 
nonlinear system may have more than one positions of static 
equilibrium. If we pick one equilibrium, x 0 we can establish 
its stability properties by linearization. The linearized 
system becomes, 

i-Ax (4 ‘ 3) 

where A is the Jacobian matrix of f (x) evaluated at x 0 , 



* ~f(x) 



(4.1) 



f(x 0 ) = 0 



(4.2) 




(4.4) 



16 



and the state has been defined to designate small deviations 
from the equilibrium x 0 , 

x — * x-x 0 (4.5) 

In system dynamics as long as all eigenvalues of A have 
negative real parts, we know that the linear system will be 
stable. This means that the equilibrium x 0 will be stable for 
the nonlinear system as well. This is in fact Lyapunov's 
linearization technique. 

B. BIFURCATION 

Values of the nonlinear system parameters at which the 
qualitative nature of the system's motion changes are known as 
critical or bifurcation values. The phenomena of bifurcation, 
i.e., quantitative change of parameters leading to qualitative 
change of system properties, is the topic of bifurcation 
theory. Euler buckling (Pitchfork bifurcation) , limit cycles 
(Hopf bifurcation) are common examples of bifurcation. 

Classical definition of stability states, that the real 
part of all the eigenvalues of the system must be negative. 
Therefore, our initial investigations into the stability of 
the SUBOFF model was to find those eigenvalues whose real 
parts cross the imaginary axis. We used the bifurcation 
analysis program, included in the Appendix, to calculate the 
eigenvalues of the system. 
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By linearizing the equations of motion, equations 2.1 
through 2.4 , the state space equations of the dynamic system 
can be written in the form, 



x =Ax+Bu 



(4.5) 



where , 



u - - K x (4.6) 

and K is the vector of controller gains, as calculated by pole 
placement in equations 3.13. The eigenvalues of the system are 
found by solving, 

det|A-BK-sI| = 0 (4.7) 

In the bifurcation program a pseudo-root locus method is 
employed where the time constant, T c , is fixed. The constant 
T c fixes to placement of the system poles at a given nominal 
forward speed U 0 and then the model speed, U, is varied 
incrementally with the system eigenvalues calculated at each 
speed increment. When the real part of an eigenvalue changes 
sign between the limits of a speed increment a bisection 
method is employed to find the speed where the real part of 
the eigenvalue equals zero. 

For each point where the real part of an eigenvalue 
crosses the imaginary axis the associated T c and U are plotted 
on a bifurcation map. This map delineates the regions of 
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classical stability (all eigenvalues on the left hand plane) 
from the regions of instability. A family of bifurcation maps 
were generated by varying nominal speed, U 0 , initial 
stability, z GB , and longitudinal center of gravity/buoyancy 
separation, x GB of the submersible. 

Figure (4.1) shows a typical bifurcation map with its five 
distinct regions [Ref. 1]. Region I is the area of classical 
stability. In region II there is one real positive eigenvalue 
which is indicative of a pitchfork bifurcation. Pitchfork 
bifurcations of this model were previously examined by Reidel 
[Ref. 1], Regions III, IV, and V have at least one pair of 
complex conjugate eigenvalues with a positive real part. This 
would indicate that there should be an unstable oscillatory 
behavior for the model. 

C. RESULTS AND DISCUSSION 

The classical stability region in bifurcation maps lies 
between pitchfork and Hopf bifurcation boundaries. The limits 
or parameters must be defined for the system designer prior to 
starting the design. By maximizing the region of stability we 
can give the designer the most leeway in his work. There are 
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TIME CONSTANT 




FIGURE 4.1 A typical bifurcation map showing the five 
distinct regions 
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two parameters that we used to change the bifurcation maps, 
the longitudinal separation between center of gravity and 
center of buoyancy, x GB and the initial stability, z GB . 

First we look at the change in x GB . In figures 4.2 through 
4.7 we plotted bifurcation curves for different initial 
stabilities as x GB varies. We can see that as x GB increases 
the Hopf bifurcation branches Hi an H2 move towards higher 
speeds and time constants and thus increasing the stability 
area. The H3 branch however remains constant. 

The other important point that we observed is that the 
system becomes unstable at nominal speed at higher time 
constants. This is unexpected because we are designing around 
our nominal speed. A more careful examination in the trimmed 
case shows that the actual forward velocity becomes, \/u 2 +w 2 . 
Therefore the system may become stable at a value of u other 
than nominal . 

The next parameter we examined was the initial stability, 
z GB . Figures 4.8 through 4.14 show the effect of varying z GB 
from .2 to .4 ft for different x GB values. The H3 branch 
remains constant while the upper speed H2 branch moves down 
effectively decreasing the area of stability. The low speed 
curve HI moves upwards and increases the low speed area of 
stability. 
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TIME CONSTANT 




FIGURE 4.2 Bifurcation map as xg changes between xg=0 and 
xg=- . 2 , zg= . 2 . 
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TIME CONSTANT 




0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 

VELOCITY (NONDIMENSIONAL) 



FIGURE 4.3 Bifurcation map as xg changes between xg=0 and 
xg=.2, zg=.2. 
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TIME CONSTANT 




0.4 0.6 0.0 1 1.2 1.4 1.6 1.8 



VELOCITY (NONDIMENSIONAL) 



FIGURE 4.4 Bifurcation map as xg changes between xg-0 and 
xg=- . 3 , zg= . 3 . 
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TIME CONSTANT 




FIGURE 4.5 Bifurcation map as xg changes between xg-0 and 
xg=.3, zg=.3. 
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TIME CONSTANT 




0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 



VELOCITY (NONDEMENSIONAL) 



FIGURE 4.6 Bifurcation map as xg changes between xg 0 and 
xg=- . 3 , zg= . 4 . 
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TIME CONSTANT 




VELOCITY (NONDIMENSIONAL) 



FIGURE 4.7 Bifurcation map as xg changes between xg=0 and 
xg=.3, zg=.4. 
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TIME CONSTANT 




FORWARD VELOCITY (NONDIMENSIONAL) 

FIGURE 4.8 The effects of changing zg on the bifurcatic 
maps , xg= . 3 . 
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TIME CONSTANT 




0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 



FORWARD VELOCITY (NONDIMENSIONAL) 

FIGURE 4,9 The effects of changing zg on the bifurcation 
maps , xg= . 2 . 
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TIME CONSTANT 




FIGURE 4.10 



The effects of changing zg on the bifurcatio 



maps, xg=.l. 
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TIME CONSTANT 




0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 

FORWARD VELOCITY (NONDIMENSIONAL) 



FIGURE 4.11 



The effects of changing zg on the bifurcation 



maps, xg=0 . 
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TIME CONSTANT 




0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 

FORWARD VELOCITY (NO ND IMENSIONAL) 



FIGURE 4.12 The effects of changing zg on the bifurcatio 
maps , xg=- . 1 . 
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TIME CONSTANT 




0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 

FORWARD VELOCITY (NONDIMENSIONAL) 



FIGURE 4.13 The effects of changing zg on the bifurcation 



maps , xg=- . 2 . 



TIME CONSTANT 




FORWARD VELOCITY (NONDIMENSIONAL) 

FIGURE 4.14 The effects of changing zg on the bifurcatio 
maps , xg=- . 3 . 
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V. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

Hopf bifurcation analysis is a very useful design tool in 
the design and evaluation phase. Hopf bifurcation analysis and 
an identification program that can evaluate the hydrodynamic 
coefficients for the submersible vehicle will be very useful 
and save money and time by reducing the amount of model 
testing. An effective set of control system parameters can be 
generated in this process that will be optimal for the final 
design of the submersible. 

This type of analysis can set the limits of the ranges of 
important parameters such as metacentric height and 
longitudinal separation of buoyancy/gravity centers. As we 
have seen changes in these two parameters can have dramatic 
effects on stability. It was found that the moderate speed 
region of stability increases with increasing metacentric 
height. The same is not true, however, for high speeds. The 
longitudinal separation of center of gravity/buoyancy can have 
a profound effect on stability. It was found that the vehicle 
may be unstable even at nominal speed. This was attributed to 
the fact that at high trim angles, the feedback gains which 
are computed at zero trim, can no longer guarantee stability. 
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B. RECOMMENDATIONS 



The bifurcation analysis program should be expanded to 
evaluate the performance of the submarine including effects of 
external forces such as wave effects, currents, and free 
surface effects. 
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APPENDIX - BIFURCATION ANALYSIS PROGRAM 



C PROGRAM BIFUR 1 . FOR 
C BIFURCATION ANALYSIS 
C PARAMETERS ARE: TC VS. U 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION Kl,K2,K3,K4,L,MQDOT,MWDOT,MQ,MW,MDS,MDB,MD, 
& MASS,IY,P1,P2,P3,P4,XGB,ZGB 

DIMENSION A(4,4),FV 1 (4),I V 1 (4),ZZZ(4,4), WR(4), WI(4),XL(25), 

& BR(25),VEC0(25),VEC1(25),VEC2(25) 

COMMON PI, P2,P3,P4 
C 

OPEN (1 1 ,FILE='B IF 1 .RES', ST ATUS =r NEW') 

OPEN ( 1 2,FILE-BIF2.RES', STATUS='NEW') 

OPEN (13, FILE-BIF3 RES', STATUS-NEW') 

C NUMERIC INFO OF DARPA SUBOFF MODEL 
WEIGHTS 556.2363 
BUO =1556.2363 
L =13.9792 
IY =561.32 
G =32.2 

MASS =WEIGHT/G 
RHO =1.94 
XB =0.0 
ZB =0.0 
CD =0.4 
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CD =0.5*CD*RHO 



C 

WRITE (*,*) 'ENTER MIN, MAX, AND INCREMENTS IN Tc (nondim)’ 
READ (*,*) TCMIN,TCMAX,ITC 

WRITE (*,*) 'ENTER MIN, MAX, AND INCREMENTS IN U (nondim)’ 
READ (*,*) UMIN,UMAX,IU 
C WRITE (*,*) 'ENTER NOMINAL SPEED' 

C READ (*,*)U0 

WRITE (*,*) 'ENTER XG AND ZG' 

READ (*,*) XG,ZG 
UO-9 
C 

ZGB=ZG-ZB 
XGB=XG-XB 
TCMIN=TCMIN*L/UO 
TCMAX=TCMAX*L/UO 
UMIN =UMIN*U0 
UMAX =UMAX*U0 
C HYDRODYNAMIC COEFFICIENTS 
ZQDOT=-6.3300E-04*0.5*RHO*L**4 
ZWDOT=-1.4529E-02*0.5*RHO*L**3 
ZQ = 7.5450E-03*0.5*RHO*L**3 
ZW =-L3910E-02*0.5*RHO*L**2 
ZDS =-5.6030E-03*0.5*RHO*L**2 
ZDB =-5.6030E-03*0.5*RHO*L**2 
MQDOT=-8.8000E-04*0.5*RHO*L !,: *5 
MWDOT=-5.6100E-04*0.5*RHO*L* ,|S 4 
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MQ =-3.7020E-03*0.5*RHO*L**4 
MW = 1.0324E-02*0.5*RHO*L**3 
MDS =-2.4090E-03*0.5*RHO*L**3 
MDB = 2.4090E-03*0.5*RHO*L**3 



XL(1)=0.0 

XL(2)=0. 1 

XL(3)=0.2 

XL(4)=0.3 

XL(5)=0.4 

XL(6)=0.5 

XL(7)=0.6 

XL(8)=0.7 

XL(9)=1.0 

XL(10)=2.0 

XL(11)=3.0 

XL(12)=4.0 

XL(13)=7.7143 

XL(14)=10.0 

XL(15)=15.1429 

XL(16)=16.0 

XL(17)=17.0 

XL(18)=18.0 

XL(19)=19.0 

XL(20)=20.0 

XL(21)=20.1 
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X L(22)=20.2 
XL(23)=20.3 
XL(24)=20.4 
XL(25)=20.4167 
DO 102 N= 1,25 

XL(N) = (L/20.)*XL(N) - L/2. 

102 CONTINUE 
BR(1)=0.0 
BR(2)=0.485 
BR(3)=0.658 
BR(4)=0.778 
BR(5)=0.871 
BR(6)=0.945 
BR(7)=1.010 
BR(8)= 1.060 
BR(9)=1.18 
BR(10)=1.41 
BR(1 1 )= 1.57 
BR(12)=1.66 
BR(13)=1.67 
BR(14)=1.67 
BR( 1 5)=1 .67 
BR(16)=1.63 
BR(17)=1.37 
BR(18)=0.919 
BR( 1 9)=0.448 
BR(20)=0. 195 



40 



BR(21)=0. 188 
BR(22)=0.168 
BR(23)=0. 132 
BR(24)=0.053 
BR(25)=0.0 
DO 104 K=l,25 
VEC0(K)=BR(K) 

VEC 1 (K)=XL(K)*BR(K) 

VEC2(K)=XL(K) * XL(K) *BR(K) 

104 CONTINUE 

CALL TRAP(25,VECO,XL,EO) 

CALL TRAP(25,VECLXL,E1) 

CALL TRAP(25,VEC2,XL,E2) 

C 

ALPHA=0.0 

ZD=ZDS+ALPHA*ZDB 

MD=MDS+ALPHA*MDB 

C CALCULATING THE SS PITCH ANGLE 0 O 
C WITH NEWTON RAPHSON METHOD 
P1=ZW*MD-MW*ZD 
P2= XG*BUO*ZD 
P3= ZG*BUO*ZD 
P4= CD*(MD*E0 - ZD*E1) 
WRITE(*,*)P1,P2,P3,P4 
EPSI=. 00000001 
THETA0=0 
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IF(XGB GT.O) P4=-1*P4 

18 DO 19 1=1,2000 
FT=FUNC(THET AO) 

DFT=DFUNC(THETAO) 

DELT=FT/DFT 

THET AO=THET AO-DELT 
IF(ABS(DELT)-EPSI) 20,20,19 

19 CONTINUE 
FT=FUNC(THET AO) 

20 WRITE(*,*) THET AO , FT 
C 

D V=(MAS S-Z WDOT) *(I Y -MQDOT)-(MAS S *XG+ZQDOT) * (MAS S * XG+MWDO T) 
A1 1DV=(IY-MQDOT)*(ZW-2*CD*EO*U*TAN(THETAO)) 

& +(MASS*XG+ZQDOT)*(MW+2*CD*E1*U*TAN(TFIETAO)) 

A12DV=(IY-MQD0T)*(MASS+ZQ+2*CD*E1 *U*TAN(THETAO))+ 

& (MASS*XG+ZQDOT) 

& *(MQ4VLASS*XG4VIASS*ZG*U*TAN(THETA0)-2*CD*E2*U*TAN(TFIETA0)) 
A1 3DV=-(MASS*XG+ZQDOT)*WEIGHT 
B 1 DV =(IY-MQDOT)*ZD+(MASS *XG+ZQDOT)*MD 
A21DV=(MASS-ZWDOT)*(MW+2*CD*E1*U*TAN(THETAO)) 

& +(MASS*XG+MWDOT)*(ZW-2*CD*EO*U*TAN(THETAO)) 

A22DV=(MASS-ZWDOT)*(MQ-MASS*XG-MASS*ZG*U*TAN(THETAO)-2*CD*E2*U 
& * T AN (THET A0))+(M AS S * XG+MWDOT) * 

& (MASS+ZQ+2*CD*E1 *U*TAN(THETAO)) 

A23 DV=-(MAS S-Z WDOT) * WEIGHT 
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B2DV=(MASS-ZWD0T)*MD+(MASS*XG+MWD0T)*ZD 

c 

A1 1=A1 1DV/DV 
A12=A12DV/DV 
A13=A13DV/DV 
A21=A21DV/DV 
A22=A22DV/DV 
A23=A23DV/DV 
B1 =B 1 D V /D V 
B2 =B2DV /DV 
C 

EPS =l.D-5 
ILMAX= 1 500 
C 

DO 1 1=1, ITC 
WRITE (*,2001) I, ITC 

TC=TCMIN+(I- 1 )*(TCMAX-TCMIN)/(ITC- 1 ) 

POLE=1.0/TC 
ALPHA3=4.0*POLE 
ALPHA2=6.0*POLE**2 
ALPHA1=4.0*POLE**3 
ALPHA0= POLE* *4 

K4=ALPHA0/((B1*A21-B2*A1 1)*U0**4+(B1*A23-B2*A13)*ZGB*U0**2) 

A2M=B1*U0**2 

A3M=B2*U0**2 

A0M=-(A1 1 + A22)*U0-ALPHA3 

B1M=B2*U0**2 
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B2M=(B2* A 1 2-B 1 * A22)*U0* *3 
B3M=(B 1 *A2 1 -B2* A1 1)*U0**3 

B0M=(A1 1*A22-A21*A12)*UO**2-A23*ZGB-AJLPHA2-B1*UO*UO*K4 
C1M=(B2*A1 1-B1*A21)*U0**3 
C2M=(B1*A23-B2*A13)*ZGB*U0**2 

C0M=(A13*A21-A23*A1 1)*ZGB*U0+ALPHA1-(B2+B1*A22-B2*A12)*K4*U0**3 

K2=C1M*B0M*A3M-B1M*C0M*A3M-C1M*B3M*A0M 

K2=K2/(C1M*B2M*A3M-B1M*C2M*A3M-C1M*B3M*A2M) 

K1=(C0M-C2M*K2)/C1M 

K3=(A0M-A2M*K2)/A3M 

DO 2 J=1,IU 

U=UMIN+(J- 1 )*(UMAX-UMIN)/(IU- 1) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(2, 1 )=-A 1 3*(XGB*SIN(THETAO)-ZGB*COS(THETAO))+B 1 *U*U*K 1 
A(2,2)=A11*U +B1 *U*U*K2 
A(2,3)=A12*U +B1*U*U*K3 
A(2,4)= B1*U*U*K4 

A(3,1)=-A23*(XGB*SIN(THETAO)-ZGB*COS(THETAO))+B2*U*U*K1 

A(3,2)= A21*U +B2*U*U*K2 

A(3,3)= A22*U +B2*U*U*K3 

A(3,4)= B2*U*U*K4 

A(4,l)=- U 

A(4,2)= 1.0D0 
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A(4,3)= O.ODO 
A(4,4)= O.ODO 
C 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 
CALL DSTABL(DEOS,WR,WI,FREQ) 

C 

IF (J.GT.l) GO TO 10 
DEOSOO= DEOS 
UOO = U 
LL= 0 
GO TO 2 

10 DEOSNN = DEOS 
UNN = U 

PR= DEOSNN*DEOSOO 
IF (PR.GT.0.D0) GO TO 3 
LL = LL+1 

IF (LL.GT.3) STOP 1000 
IL = 0 
UO = UOO 
UN = UNN 
DEOSO = DEOSOO 
DEOSN = DEOSNN 
6 UL = UO 
UR = UN 
DEOSL = DEOSO 
DEOSR = DEOSN 
C U = (UL+UR)/2.D0 
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ALPHA = (DEOSL-DEOSR)/(UL-UR) 

U =- (DEOSL-ALPHA*UL)/ ALPHA 

A(l,l) = O.ODO 

A(l,2) = O.ODO 

A(l,3) = 1.0D0 

A(l,4) = O.ODO 

A(2,l) = A13*(XGB*SIN(THETAO)-ZGB*COS(THETAO))+Bl*U*U*Kl 
A(2,2) = A1 1*U +B1*U*U*K2 
A(2,3 )= A12*U +B1 *U*U*K3 
A(2,4 )= B1 *U*U*K4 

A(3,l) = A23 * (XGB * SIN(THET AO)-ZGB * COS(THET AO))+B2* U* U*K 1 

A(3,2) = A21*U +B2*U*U*K2 

A(3,3) = A22*U +B2*U*U*K3 

A(3,4 ) = B2*U*U*K4 

A(4, 1 ) =- U 

A(4,2) = l.ODO 

A(4,3) = O.ODO 

A(4,4 )= O.ODO 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL(DEOS,WR,WI,FREQ) 

DEOSM = DEOS 
UM = U 

PRL = DEOSL*DEOSM 
PRR = DEOSR*DEOSM 
IF (PRL.GT.O.DO) GO TO 5 
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UO = UL 
UN = UM 
DEOSO = DEOSL 
DEOSN = DEOSM 
IL = IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF = DAB S(UL-UM) 

IF (DIF.GT.EPS) GO TO 6 
U = UM 
GO TO 4 

5 IF (PRR.GT.0.D0) STOP 3200 
UO = UM 
UN = UR 
DEOSO = DEOSM 
DEOSN = DEOSR 
IL = IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF = DABS(UM-UR) 

IF (DIF.GT.EPS) GO TO 6 
U = UM 

4 LLL = 10+LL 

WRITE (LLL,*) U/U0,TC*U0/L 
3 UOO = UNN 

DEOSOO = DEOSNN 
2 CONTINUE 
1 CONTINUE 
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2001 FORMAT (215) 
END 



SUBROUTINE DSTABL(DEOS,WR,WI,OMEGA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION WR(4),WI(4) 

DEOS=-1.0D+20 
DO 1 1= 1,4 

IF (WR(I).LT.DEOS) GO TO 1 
DEOS = WR(I) 

IJ = I 

1 CONTINUE 
OMEGA = WI(IJ) 

OMEGA = DABS(OMEGA) 

RETURN 

END 

C 

SUBROUTINE TRAP(N,A,B,OUT) 

C NUMERICAL INTEGRATION ROUTINE USING THE TRAPEZOIDAL RULE 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION A(1),B(1) 

N1 = N-l 
OUT = 0.0 
DO 1 1= 1,N1 

OUT1 = 0. 5 *(A(I)+A(I+ 1 ))*(B(I+ 1 )-B(I» 

OUT = OUT+OUT1 
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I CONTINUE 
RETURN 
END 

C FUNCTIONS USED IN NEWTON RAPHSON ROUTINE 
FUNCTION FUNC(THETAO) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON PI, P2,P3,P4 

FUNC = P1*DTAN(THETA0) + P2*DCOS(THETAO) + P3*DSIN(THETA0) 

& + P4*(DTAN(THETA0))* t 2 

RETURN 
END 

FUNCTION DFUNC(THETAO) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON PI, P2,P3,P4 

DFUNC= P1*(LDO/DCOS(THETAO))**2 - P2*DSIN(THETA0) + 

& P3*DCOS(THETAO) + P4*2.DO*DTAN(THETAO)*(1.DO/DCOS(THETAO))**2.DO 
RETURN 
END 
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